#ifndef UHM_DISCRETE_EXTENSION_METIS_HXX_
#define UHM_DISCRETE_EXTENSION_METIS_HXX_


namespace uhm {


  // Metis graph format
  template<class M_>
  Bool_ elements2graph(std::vector< typename M_::Element_* > *els,
                       std::vector<Int_> *xadj,
                       std::vector<Int_> *adjncy,
                       std::vector<Int_> *adjwgt,
                       Int_ pid=0) {

    Int_ s = (Int_)els->size();

    time_in("uhm:: elements2graph:: Creating Adjancy Matrices",UHM_TIME_LEVEL_INTERNAL);

    // create subset index and mark its parent
    for (Int_ i=0;i<s;++i) 
      els->at(i)->set_marker(i);

    // create graph
    std::vector< std::map< Int_, Int_ > > adj; 
    adj.assign(s, std::map<Int_,Int_>());

    // for all elements
#pragma omp parallel for schedule(static)
    for (Int_ i=0;i<s;++i) {
      typename M_::Element_* e = els->at(i);

      // loop through nodes
      typename std::vector< typename M_::Element_::Entry_ >
        ::iterator nit;
      for (nit=e->nds().begin();nit!=e->nds().end();++nit) {

        // exclude boundary nodes
        typename M_::Node_* n = nit->node();        
        if (n->get_kind() == M_::NODE_DEFAULT) {

          // loop owners
          typename std::set< typename M_::Element_* >::iterator eit;
          for (eit=n->owner().begin();eit!= n->owner().end();++eit) {

            // if owner has partition info matching pid, add it in adjacency
            typename M_::Element_* o = (*eit);
            if (o->id().get_part()==pid && o!=e) {
              
              // vertex point and its weight
              std::pair< std::map< Int_, Int_ >::iterator, Bool_ > ret;
              std::pair< Int_,Int_ > in;
              in.first  = o->get_marker();
              in.second = (n->id().get_dof() + n->id().get_weight());

              ret = adj[i].insert(in);

              if (ret.second==false)
                ret.first->second += in.second;
            }
          }
        }
      }
    }
    time_out(UHM_TIME_LEVEL_INTERNAL);

    // clear input array
    if (xadj->size()) xadj->clear();
    if (adjncy->size()) adjncy->clear();
    if (adjwgt->size()) adjwgt->clear();


    time_in("uhm:: elements2graph:: Converting Compressed Format",UHM_TIME_LEVEL_INTERNAL);

    xadj->reserve   ( s+1 );
    xadj->push_back ( 0 );
    for (Int_ i=0;i<s;++i) {

      std::map<Int_,Int_>::iterator ait;
      for (ait=adj[i].begin();ait!=adj[i].end();++ait) {
        adjncy->push_back(ait->first);
        adjwgt->push_back(ait->second);
      }
      xadj->push_back(adjncy->size());
    }

    time_out(UHM_TIME_LEVEL_INTERNAL);

    /*
    std::cout << " Priting graph\n";
    for (int i=0;i<xadj->size();++i)
      std::cout << xadj->at(i) << " ";
    std::cout << std::endl;

    for (int i=0;i<adjncy->size();++i)
      std::cout << adjncy->at(i) << " ";
    std::cout << std::endl;

    for (int i=0;i<adjwgt->size();++i)
      std::cout << adjwgt->at(i) << " ";
    std::cout << std::endl;
    */

    return true;
  }

  
  inline std::ostream& graph_show(std::ostream& os,
                                  std::vector<Int_> &xadj,
                                  std::vector<Int_> &adjncy,
                                  std::vector<Int_> &adjwgt) {
    for (Int_ i=0;i<(xadj.size()-1);++i) {
      os << "************* xadj " << xadj[i] << " -- " << xadj[i+1] << std::endl; 
      for (Int_ j=0;j<(xadj[i+1] - xadj[i]);++j) 
        os << adjncy[xadj[i]+j] << " -- " << adjwgt[xadj[i]+j] << std::endl;
    }
    return os;
  }


}

#endif
